The insulator/Chern-insulator transition in the Haldane model 
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We study the behavior of several physical properties of the Haldane model as the system un- 
dergoes its transition from the normal-insulator to the Chern-insulator phase. We find that the 
density matrix has exponential decay in both insulating phases, while having a power-law decay, 
more characteristic of a metallic system, precisely at the phase boundary. The total spread of the 
maximally-localized Wannier functions is found to diverge in the Chern-insulator phase. However, 
its gauge-invariant part, related to the localization length of Resta and Sorella, is finite in both 
insulating phases and diverges as the phase boundary is approached. We also clarify how the usual 
algorithms for constructing Wannier functions break down as one crosses into the Chern-insulator 
region of the phase diagram. 



PACS numbers: 73.43.-f, 73.20.At, ll.30.Rd 



I. INTRODUCTION 

The band structure of any insulator is characterized by 
a certain discrete topological index known as the Chern 
invariant^ which encodes information about the phase 
evolution of the Bloch functions around the boundary of 
the Brillouin zone (see Sec. II). Insulators can thus be 
classified as "normal insulators" or "Chern insulators" 
depending on whether or not the Chern invariant van- 
ishes. The latter case requires breaking of time-reversal 
symmetry, so insulating ferromagnets and ferrimagnets 
could be candidates for Chern insulators. While models 
for Chern insulators can be constructed theoretically,'^ no 
experimental realizations are yet known to occur in na- 
ture. A Chern insulator, if found to exist, would have 
the remarkable feature of showing a quantum Hall effect 
in the absence of a macroscopic magnetic field. Hence, 
Chern insulators may also be referred to as "quantum 
Hall insulators." 

Although the basic theory of Chern insulators was for- 
mulated in the 1980's, not much is known theoretically 
about the general features of the electronic band struc- 
ture of such insulators. In the last 15 years or so, the 
theory of normal insulators has been greatly enriched by 
a deeper understanding of electric polarization,^ orbital 
magnetization,*'^ linear-scaling theory, methods for con- 
structing Wannier functions,^ the spatial decay of Wan- 
nier functions and of the one-particle density matrix,^ 
and related measures of localization.®"^" However, all of 
this work implicitly assumed the presence of time-reversal 
symmetry, and thus was limited to the case of normal in- 
sulators. It is therefore of considerable interest to revisit 
many of these same issues, and to reconsider whether, 
or how, the previous conclusions generalize to the case 
of Chern insulators. For example, what are the decay 
properties of the one-particle density matrix in a Chern 
insulator? Can Wannier functions be constructed, and if 
not, in what way do the usual construction procedures 
fail? If one inspects closely related measures of local- 
ization such as the gauge-invariant part of the Wannier 



spread functional,^ the localization length of Resta and 
Sorella,* or the second-cumulant moment of the electron 
distribution,^' does the localization remain finite in a 
Chern insulator, or does it diverge? 

Furthermore, an intriguing feature of Chcrn-insulating 
systems is that a phase boundary separating the Chern 
insulator from a normal insulator may occur. Such a 
normal-insulator/Chern-insulator (NI/CI) transition is 
an example of a class of topological transitions that have 
become of considerable current interest, at which a topo- 
logical invariant changes discontinuously across a phase 
boundary. Such transitions normally appear within the 
theory of correlated states. ^^'^"^ The NI/CI transition, on 
the other hand, occurs in a non-interacting context, and 
can therefore be studied at a level of detail, and tested 
with numerical calculations, in a way that is difficult 
for correlated models. In addition, Chern insulators are 
closely related to so-called "spin Hall insulators," which 
have also been the subject of a recent surge of interest.^* 
Thus, there is a special opportunity associated with the 
study of this particular topological insulator/insulator 
transition. 

In this paper, we investigate several aspects of the elec- 
tronic structure near the NI/CI transition in the two- 
dimensional Haldane model. ^ We choose the Haldane 
model because it is one of the simplest models that ex- 
hibits a quantum-Hall insulator state. The underlying 
idea of this model is to break time-reversal symmetry so 
that the transverse conductivity axy, which is odd under 
time-reversal, can become nonzero. Usually, the quan- 
tum Hall effect is associated with a gap at the Fermi 
level resulting from a splitting of the spectrum into Lan- 
dau levels by a macroscopic magnetic field. In the Hal- 
dane model, however, there is a degeneracy between the 
valence and conduction bands at certain high-symmetry 
k-points when both inversion and time-reversal symme- 
try are present. If a gap is opened by the breaking of 
inversion symmetry, the system becomes a normal insu- 
lator. However, if the gap opens as a result of breaking 
time-reversal symmetry, the system turns into a Chern 
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insulator. 

We have organized this paper as foUows. In Sec. II we 
introduce the Chern invariant, which will henceforth be 
used to classiiy the state of our system. The basics of 
the Haldane model are reviewed in Sec. III. Thereafter, 
we focus on the problems occurring when constructing 
Wannier functions for Chern insulators (Sec. IV), the be- 
havior of the spread functional (Sec. V), and the decay 
of the density matrix (Sec. VI) as the system transitions 
from the normal insulator phase into the Chern-insulator 
phase. We conclude and give an outlook in Sec. VII. 

II. THE CHERN INVARIANT 

We restrict ourselves to the case of a one-particle 
Hamiltonian H having Bloch eigenvalues and eigen- 
states IV'nk)- The cell-periodic part of the Bloch function 
Mnk;(r) = e~*'' ''?/^„k(r) is then an eigenfunction of the ef- 
fective Hamiltonian i?(k) = e~^^''^ H e^^''^ . We consider 
electrons to be spinless, but factors of two can easily be 
inserted for non-interacting spin channels. 

We can now define the Chern invariant^ for an insu- 
lator, defined here as a system with a gap in the single- 
particle density of states separating occupied and unoc- 
cupied states, to be 

occ 

C = — dk y^idkUnul X |5ku„k) , (1) 

where BZ denotes an integral over the Brillouin zone 
and dk = d/dk. The cross product notation in Eq. (1) 
implies, for example, that Cz contains terms involving 
{dk^UnuldkyUn-k) - {dkyU„\^\dk^Un\c) ■ For non- interacting 
electrons, ^^'^^ the Chern invariant is quantized in units 
of reciprocal-lattice vectors G. For the case of a two- 
dimensional system with only a single occupied band, 
Eq. (1) becomes 

C=^ [dk (^kUkl X l^kWk) . (2) 

In two dimensions the Chern invariant is a pseudo-scalar 
called the Chern number which can only take integer val- 
ues. Alternatively, we can write the Chern number in 
terms of the Berry connection A(k) = i(Mk|i9k|uk) and 
the Berry curvature f](k) = Vk x A(k) as 

C = -!- [ dk n{k) = ^ idk- A(k) . (3) 

A Chern insulator is now simply defined as an insulator 
with a nonzero Chern invariant. Conversely, we define 
a normal insulator to be an insulator with zero Chern 

invariant. Hence, the NI/CI transition is characterized 
by a change of the Chern invariant from zero to a nonzero 
value. 

The Chern invariant of Eqs. (1) and (2) is gauge 
invariant,^ i.e., invariant with respect to the choice of 




FIG. 1: (a) Four unit cells of the Haldane model. Filled 
(open) circles have site energy —A (-I-A). The first-nearest 
neighbor hopping ti is real, while the second-nearest neigh- 
bor hopping t2e"'' has a complex phase. Arrows indicate the 
direction of positive phase hopping. The Wigncr-Scitz unit 
cell corresponds to the hexagon in the center of the plot, (b) 
First Brillouin zone of the Haldane model with high-symmetry 
points marked. 



phases of the \urik), or in the more general multiband 
case, the choice of unitary rotations applied to transform 

the occupied states among themselves at a given k. It can 
be shown that in normal insulators it is always possible 
to make a gauge choice such that the Bloch orbitals are 
periodic in k-spacc (i.e., IV^nk+c) = IV-'nk)) and smooth 
in k (i.e., continuous and differentiable), whereas no such 
choice is possible for a Chern insulator.^'' 



III. THE HALDANE MODEL 



Here we provide a brief review of Haldane's model and 
its properties, as discussed in detail in Ref. [2]. As il- 
lustrated in Fig. 1. the Haldane model is comprised of 
a honeycomb lattice having two tight-binding sites per 
cell with site energies ±A, a real first-neighbor hopping 
ti, and a complex second-neighbor hopping t2e^*'^. The 
model can also be thought of as consisting of two sub- 
lattices A and B corresponding to the sites with energies 
-I- A and —A, respectively. Note that the macroscopic 
magnetic flux through the unit cell is indeed zero, result- 
ing in a vanishing macroscopic magnetic field. This fol- 
lows directly from the fact that the first-nearest neighbor 
hopping is real and no phase is picked up when hopping 
aroimd the Wigncr-Scitz imit cell. This, however, does 
not rule out a microscopic magnetic field that averages 
to zero over the unit cell. Note that the wavevector k is 
still a good quantum number under these conditions. 

Let ai,a2, and as be the vectors pointing from a site 
of the B sublattice to its three nearest A neighbors, such 
that z • ai X a2 > and x • ai > 0. If we furthermore 
define the vectors bi = a2 — aa, b2 = aa — ai, and 
ba = ai — a2, then the Hamiltonian of the Haldane model 
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(p [in units of TC] marked in Fig 2. The inset shows a magnification of the bands 

at K. Note that at (A/t2)cr the dispersion is linear. 



FIG. 2: Chern number of the bottom band of the Haldane 
model as a function of the parameters ip and A./t2 {ti = 1, 
t2 = 1/3). Analytic results are plotted as solid sinusoidal 
lines, whereas numerical results are depicted by circles. The 
vertical line shows the range of parameters that we have cho- 
sen for all our calculations. 

can be written as 

if(k) = I 2t2 cos y cos(k • hi) 

i 

+ CTi ti cos(k • SLi) + (72 ti sin(k • a^) 

i i 

+ <73 (a - 2t2 sin ip sin(k • bj)) , (4) 

i 

where the cr^ are the Pauli matrices and I is the identity 
matrix. 

The Chcrn number can now be calculated analytically 
or numerically according to Eq. (2) or (3). For our tests, 
we have chosen a lattice constant equal to unity, ti = 1, 
and f2 = 1/3. If the Chern number of the bottom band is 
mapped out as a function of the remaining model param- 
eters If and A/t2, wc obtain the Haldane phase diagram 
shown in Fig. 2. Since we are interested in studying the 
transition from a normal insulator to a Chern insulator, 
we choose for all our calculations below a path in the 
phase diagram that crosses the phase boundary. Specif- 
ically, we traverse the vertical line in Fig. 2 where the 
phase (f is fixed at 7r/4 and A/t2 is reduced from 6 to 2. 
At the critical value (A/i2)cr — 3\/3 sin(7r/4) w 3.67, the 
phase boundary is crossed. 

The band structure of the Haldane model is plotted in 
Fig. 3 along some high-symmetry lines in the Brillouin 
zone (see Fig. lb). It shows a remarkable feature as the 
system passes through (A/i2)cr- In the normal-insulator 
region, the two bands are separated by a finite gap. As 
the critical value is approached, the gap at K gets smaller 
and smaller. Finally, exactly at (A/i2)cr the bands touch 
at K in such a way that the dispersion relation is linear. 
Such points are also referred to as Dirac points. When 
going further into the Chern-insulator region, the bands 



separate again. Note that our specific choice of ti = 1 
and t2 = 1/3 prevents the bands from overlapping. If A 
and t2 sin ip are both chosen to be zero, two Dirac points 
form at K and K', and the Haldane model then becomes 
an appropriate model for a graphene sheet. 

In the normal-insulator region of the Haldane model 
the Chern number of each band is zero, so that the to- 
tal Chern number (the sum of the Chern numbers of the 
upper and lower bands) is obviously also zero. When 
the phase boundary is crossed, the Chern numbers of the 
upper and lower bands become ±1, but their sum still 
remains zero. The closure and reopening of the gap as 
the NI/CI boundary is crossed corresponds to the "dona- 
tion" of a Chern unit from one band to another through 
the temporarily formed Dirac point. In the present case, 
the total Chern number must always remain zero because 
the model, having a tight-binding form, assumes Wannier 
representability of the overall band space, and a non-zero 
Chern number is inconsistent with such an assumption. 
More generally, the total Chern number of a group of 
bands should not change when a gap closure and reopen- 
ing occurs among the bands of the group, as long as the 
gaps between this group and any lower or higher bands 
remains open. 

It is possible to argue on very general grounds that a 
finite sample cut from a Chern insulator must have con- 
ductive channels, otherwise known as chiral edge states, 
that circulate around the perimeter of the sample^^ in 
much the same way as for the quantum Hall effect. ^"'^^ It 
is therefore of interest to investigate the electronic struc- 
ture of the Haldane model from the point of view of the 
surface band structure. We consider a sample that is fi- 
nite in the ba direction (specifically, 30 cells wide) and 
has periodic boundary conditions along the b2 direction 
(the hi are defined above Eq. (4)); its states can be la- 
beled by a wavevector ky running from —w/a to -|-7r/a, 
where a is the repeat unit in the y direction. The en- 
ergy eigenvalues are plotted vs. ky for several values of 
A/t2 in Fig. 4. At first sight, the surface band struc- 
ture shows qualitatively the same information as the bulk 
band structure in Fig. 3. For A/t2 = 6, the valence and 
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FIG. 4: Energy vs. wavevector ky for the Haldane model in a 
strip geometry 30 cells wide along ba direction and extending 
infinitely along b2 direction. For A/t2 < (A/t2)cr (bottom 
panel), chiral edge states are visible. 

conduction bands are separated by a finite gap. At the 
Chern transition a Dirac point forms, showing the char- 
acteristic linear dispersion expected around such a point. 
Howcvcir. when we go deeper into the Chern insulator, 
the surface band structure reveals a new behavior: one 
surface band now crosses from the lower manifold to the 
upper one with increasing ky^ and another crosses in the 
opposite direction. Further inspection shows that the up- 
going and downgoing states are localized to the right and 
left surfaces of the strip, respectively. Thus, if the Fermi 
level lies in the bulk gap, there will bo metallic states with 
Fermi velocities parallel to the surfaces and with oppo- 
site orientation, i.e., a chiral (counterclockwise) circula- 
tion of edge states around the perimeter of the sample, 
as expected. 



IV. BREAKDOWN OF WANNIER-FUNCTION 
CONSTRUCTION AT THE CHERN 
TRANSITION 

We now study aspects of the NI/CI transition that are 
related to Wannier functions (WFs) and electron local- 
ization. We expect that in the normal-insulator phase, it 
should be straightforward to construct Wannier functions 
via a k-space construction. The term "Wannier function" 
is usually applied only in the case of periodic systems, but 
for finite samples one can construct well localized Boys 



orbitals^^ which play the same role and which map onto 
the WFs in the thermodynamic limit n — > oo. Thus, if 
wo cut a finite sample from a normal-insulator realization 
of the Haldane model, we also expect it to be straightfor- 
ward to construct such Boys orbitals. The question then 
arises as to what, precisely, will "go wrong" with these 
procedures if one tries to do the same on the Chern- 
insulator side of the transition. In particular, for a finite 
sample cut from the Haldane model, it is unclear how 
the system would "know" whether the finite sample cor- 
responds to the normal-insulator or Chern-insulator side 
of the transition, and how the construction would break 
down in the latter case. In this Section, we investigate 
these issues, first in the context of the real-space con- 
struction, and then later from the k-space point of view. 

We start, then, by considering finite n x n samples of 
the Haldane model. We can interpret Fig. 1 as show- 
ing a picture of a finite sample of size n = 2; we study 
similarly-constructed samples of size n =10, 20, 30, and 
40. For each sample, the Boys orbitals are constructed 
as follows. We define the projection operator onto the 
occupied states as 

occ 

P = Y,\i>n){i>n\, (5) 

n 

and we choose a set of well- localized "trial" orbitals \ta), 
equal in number to the number of occupied states, that 
we want the Boys orbitals to be roughly modeled af- 
ter. We then construct the projected trial functions 
Ij/q) = P\ta)- Since p(r,r') = (r|P|r') is expected 
to decay exponentially in |r — r'| for an insulator (see 
Sec. VI), wc expect the \ya) to be localized as well, and 
as long as they are not overcomplete they will span the 
occupied space of interest. However, they are not or- 
thonormal, so the last step is to carry out a symmet- 
ric orthonormalization.^^ This is done by computing the 
overlap matrix 

Sap = {Valyp) (6) 

and then constructing the final Boys orbitals \LOa) as 
h„)=^(5-V2)^„|y^). 

In the context of the Haldane model, it is natural to 
choose the trial functions to be a set of (5-functions lo- 
cated on the sites of the lower-energy sublattice. With 
this choice, we can now study the lowest and highest 
eigenvalue of Sap as the parameter A/t2 traverses the 
path shown in Fig. 2. While the highest eigenvalue 
remains very close to 1, the lowest eigenvalue drops 
and rapidly approaches zero in the Chern-insulator re- 
gion, i.e., for A/t2 values below the critical value of 
(A/t2)cr ~ 3.67, as shown in Fig. 5. The slope of the 
"drop" depends on the size of the sample and becomes 
steeper as the sample size gets larger. For any given 
value of A/<2 < (A/t2)cr, the lowest eigenvalue appears 
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FIG. 6; Plot of s(k) of Eq. (12) along some high-symmetry 
lines for several values of A/t2. 



FIG. 5: Lowest eigenvalue of the overlap matrix Safi = 
{ya\yi3) as a function of the parameter A/t2 for different sam- 
ple sizes n X n. Note the logarithmic scale. 

to approach zero exponentially with sample size. When 

the eigenvalue becomes too small, the inversion to ob- 
tain S^^^^ becomes ill-conditioned, and the symmetric 
orthonormalization in Eq. (7) can no longer be carried 
out. It follows that Boys orbitals cannot be constructed 
in the Chern-insulator phase, at least not using this ap- 
proach. 

We now change perspective and look at the problem 
from the k-space point of view, where we find that some- 
thing similar happens. WFs for periodic samples are de- 
fined by 

1^") = 7^ / '^^ e-*-«-|V„k) , (8) 

\^'^) Jbz 

where the inverse relation is 

|V'nk)=^e^'^-^|Rn). (9) 

R 

In this notation |Rn) refers to the n'th WF in cell R. 

As mentioned previously, for systems with zero Chern 
invariant, the Bloch orbitals can always be chosen to obey 
a smooth and periodic gauge IV'nk+c) = IV'nk). However, 
if the Chern invariant becomes nonzero, this choice is no 
longer possible. In this case it is possible to make a 
periodic gauge choice that is smooth almost everywhere, 
but there must be singularities ( "vortices" ) somewhere in 
the interior of the BZ. For example, in two dimensions, 
assume a gauge choice that is periodic and also smoothly 
defined everywhere in the BZ except in a small disk lo- 
cated somewhere in the interior of the BZ. The periodic 
gauge choice implies that / dk ■ A(k) around the perime- 
ter of the BZ must vanish. Applying Stoke's theorem as 
in Eq. (3), but now to the region excluding the small disk, 
implies that § dk ■ A(k) around the circumference of the 
small disk must approach — C in the limit that the disk 
becomes small. For the Chern phase (C 7^ 0), this im- 
plies that there must be a vortex singularity in the phase 
choice inside the disk. If one attempts to construct WFs 
naively using Eq. (8), one then finds that the discontinu- 
ity in the phase choice of |unk) at the vortex in k-space 



leads to the destruction of exponential localization of the 
WFs in real space. 

We have investigated how this problem manifests it- 
self if one attempts to construct WFs using standard 
k-space methods. Similar to the approach described in 
Ref. [6], we again adopt a projection method in which one 
chooses trial Bloch-like functions |tk) that are smooth 
and periodic in k-space. This can be done by construct- 
ing the \tk) from a set of real-space trial functions \ta), 
i.e., t]f,{r) = X^R,e''*'^ta(r — R). Then one can construct 
projected states Ij/k) via 

\yu) = P\t^.) = J2 IMi^M = IV'k)(Vk|ik) (10) 

k' 

and orthonormalized projected states 
where 

s(k) = (yk|2/k) = |(tk|yk)|'. (12) 

The WFs are then constructed by Fourier transforming 
to real space using Eq. (8) with \w\s) substituted for iV'k)- 
Clearly, if s(k) should vanish at some k, this procedure 
would fail. 

We can now study what happens if this construction 
procedure is applied to the Haldane model. We again use 
trial functions that are ^-functions located on the lower- 
energy sites. We study the behavior of s(k) as a function 
of k throughout the BZ, while varying A/t2 along the 
line in Fig. 2. Results for s(k) are plotted along some 
high-symmetry lines in Fig. 6. In the normal-insulator 
region, we find < s(k) < 1 for all k. After the phase 
boundary has been crossed at (A/t2)cr7 we find that there 
is one point k^ in the BZ for which s{ka) = 0. There is 
also one point kf, for which s(kf,) = 1 exactly. In our nu- 
merical calculations, the locations of ka and k;, coincide 
with the points K and K' respectively. By experiment- 
ing with different trial functions, we have found that the 
precise locations of the minimum and maximum may de- 
viate from K and K', and the value at the maximum may 
be less than unity. However, we always find a point k^ at 
which s{ka) — 0. This is the point at which ('(/'kl^k) = 0; 
the robustness of such a zero-crossing can be understood 
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FIG. 7: Density of s(k) values as a function of A/t2 obtained 
for periodic samples with a 300 x 300 k-mesh. 



heuristically by realizing that by adjusting the two pa- 
rameters kx and ky, the real and imaginary parts of the 
complex scalar (^/'kl^k) can gcnerically both be made to 
vanish. From Eqs. (10 12) it follows that the phase of 
|t/k) evolves by 27r as one circles around ko, so that a 
vortex-like singularity is generated in the phase of |w;k) 
about ka, with |wk) becoming ill-dcfincd precisely at ka. 
Thus, the construction of well-localized WFs is no longer 
possible. 

Instead of focusing only on the lowest eigenvalue, we 
plot in Fig. 7 the "density of overlap values" s(k). In 
the normal-insulator region of A/t2, one sees typical 2D 
van Hove singularities, and in particular, a well-defined 
minimum above zero. In the Chern-insulator region of 
A/t2, on the other hand, the density of overlap values 
shows a tail extending all the way to zero. 

In summary, when the system is in its normal-insulator 
phase, the construction of Boys orbitals for finite sam- 
ples, or of WFs for periodic samples, can be carried out 
in the usual way using a projection method. However, 
once the NI/CI phase boundary has been crossed, such 
a construction is bound to fail because of singularities 
that appear in the overlap matrices in both the real-space 
finite-sample and k-space extended-sample approaches. 



V. THE SPREAD FUNCTIONAL 

Another quantity that shows interesting behavior as 
the phase boundary is crossed is the spread functional ft 
in real space, defined by Marzari and Vanderbilt (MV)^ 
to be 

O = ^ [(On|r^|On) - (0n|r|0n)2 , (13) 



where |0n) refers to the WF |Rn) for band n in the home 
unit cell R = and the sum is over occupied bands 
of the insulator. The spread functional is a measure of 
how "spread out" or delocalized the WFs are. In the 
remainder of this section, we specialize for simplicity to 



the case of a single band in two dimensions, so that fl = 
(0 1 1 0) - (0 1 r 1 0) 2 . MV showed that the spread functional 
can be decomposed as f2 = il/ -|- O, where 



O, = (0|r2|0) 



E 

R 



|(0|r|R) 



and 



9. 



(14) 



(15) 



are gauge-invariant and gauge-dependent contributions, 
respectively. The gauge-invariant part has been shown 
to be a useful measure for characterizing the system: f2/ 
is finite in insulators and diverges in metals.*^ 

MV also gave corresponding k-space expressions for 
the two parts of the functional. Defining the metric ten- 



sor : 
(and 5^ 
as 



and 



Re (9^Uk|Qk|9yUk) where Qk 



1 



|Mk)(uk| 



d/dk,j), these two quantities can be rewritten 



A 



dkTr[<?(k)- 

BZ 



A 



(2;r)2 



dk |A(k) 

BZ 



(16) 



(17) 



where A is the unit cell area, Tr [g\ = gxx + Syy, and A 
is the BZ average of A(k) defined just above Eq. (3). 

In the case of a Chern insulator, the use of the real- 
space expressions (14-15) becomes problematic, since 
well-localized WFs cannot be constructed. Nevertheless, 
the reciprocal-space expressions (16-17) remain well- 
defined. It is interesting, then, to see how these quantities 
behave in a Chern insulator. Do each of these quantities 
remain finite, or does one or both of them diverge? Also, 
what is the behavior of these quantities as one approaches 
the NI/CI phase boundary? 

To answer these questions, we have computed the 
quantities in Eqs. (16-17) using the finite-difference ver- 
sions of these equations given in Eqs. (34) and (36) of 
Ref. [6] . For the calculation of the gauge-dependent part 
VL, we have fixed our gauge such that ji/'k) is real for all 
k on the lower-energy site in the home unit cell. The 
results are plotted in Fig. 8 for different densities of the 
k-mesh. It can be seen that VLj is finite inside both the 
normal and Chern-insulator regions. At the critical value 
of (A/t2)cr ~ 3.67, however, fi/ diverges logarithmically 
with the number of k-points. Furthermore, f2 is finite in 
the normal insulator region, but diverges logarithmically 
with the number of k-points for Chern insulators. This 
latter behavior is consistent with the presence of a vortex 
in the phases of the |wk) aromid point ka, which causes 
A to diverge as |k — ka|~^ and imparts a logarithmic 
divergence to Eq. (17). It follows that the total spread 
ri is finite in normal insulators and divergent in Chern 
insulators. Heuristically, it is tempting to associate this 
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FIG. 8: Gauge-independent part flj and gauge-dependent 
part f2 of the spread functional for the Haldane model as 
a function of the k-mesh density. 



divergence with the presence of the metallic chiral edge 
states that are required to exist in Chern insulators (see 
Sec. Ill), but it is unclear precisely how these features are 
related. Note that electron localization in the quantum 
Hall regime is discussed in detail by R. Resta in Ref. [24] . 



VI. DECAY OF THE DENSITY MATRIX 

The decay of the density matrix is a fundamental prop- 
erty of a system and it is closely connected to the elec- 
tron localization. It was first studied by W. Kohn for 
one-dimensional insulators, and many others have in- 
vestigated this topic thereafter. 

For periodic samples the density matrix is defined as 



P(r,r') 



A ^ f 

— Y^ dk Ck(r)V„k(r') 



(2' 



(18) 



where we assume that the wave functions ipnk are nor- 
malized to one unit cell of area A. If the wave functions 
are written in terms of some basis functions (/'^(r), 



Vnk(r) = 5]C^„0^(r), 



(19) 



this becomes 



^ ' n=l a/3 -^^2 

(20) 

The C^^ arc the eigenvectors obtained by diagonalizing 
the model Hamiltonian — in our case Eq. (4). In a tight- 
binding model, the basis functions <?ia(r) are made up of 
localized orbitals cj) at sites Ta'. 

W = E - (R + ra)) • (21) 

R 



ax 

a 

^ -20 



-30 






1 , 1 , 

^11 











^12 




20 40 60 80 

R=Rx [lattice const.] 



40 60 80 

R = X [lattice const.] 



FIG. 9: Logarithm of the density-matrix kernel \^ai3{R^)\ of 
the Haldane model for several values of A/t2- The linear 
asymptotic behavior indicates exponential decay except at 
(A/t2)cr 3.67 (solid curves) where the decay is power-law 
instead. 



Inserting Eq. (21) into Eq. (20) gives 

p(r,r') = J2 ^a/3(R'-R)</'*(r-R-rJ</.(r'-R'-r^) 

(22) 



RR' 

a/3 



where 



A f 
(27r) ^Jbz 

The density matrix cannot be evaluated explicitly with- 
out the knowledge of the orbitals 0, but we can study 
instead the decay of ^a^CR), which essentially has the 
interpretation of being a density matrix expressed in a 
tight-binding representation. 

Calculating the decay of ^q,^(R) in Eq. (23) numer- 
ically is very demanding and the corresponding results 
are to be interpreted with caution. To ensure high accu- 
racy, we used a very dense k-mesh of 2000 x 2000 points 
and 128-bit arithmetic. Results for ^a/3(-Rx) (i.e., along 
the X direction) for the Haldane model are collected in 
Fig. 9. In normal insulators the density matrix decays 
exponentially with a power-law prefactor.^ We therefore 
choose to fit our results according to ^ai3 ~ R~°'e~^^, 
where i? = |R| = |i?x|, and a and b are fit parame- 
ters. More specifically, we performed least-square fits of 
ln|^„^| for distances up to 100 unit cells. For the decay 
behavior at the NI/CI boundary, we even went as far as 
500 unit cells. 

Within fitting error, the best-fit values for the param- 
eter b are the same for all ^ai3- Numerical results cor- 
responding to A/t2 values of 6, 5, 4, 3.67, 3, and 2 
are 0.69±1, 0.43±1, 0.118±5, O.OOOlil, 0.282±1, and 
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0.75±1, respectively. In general, when approaching the 
phase boundary from either side, the best-fit value of 
the parameter b decreases and takes its minimum of zero 
at (A/f2)cr- In other words, in the normal and Chern- 
insulator regions the decay is dominated by the exponen- 
tial behavior. However, exactly at the phase boundary 
the exponential decay vanishes {b = 0) and a pure power- 
law behavior remains, similar to metals. At (A/t2)cr the 
power-law decay is then characterized by a = 3.01±3 for 
^11 and ^22i and a = 2.00±2 for and ^211 which sug- 
gests that the "true" values are the integers 3 and 2. Note 
that the results depicted in Fig. 9 correspond to a partic- 
ular direction in real space (R — Rx.). While the decay 
parameters inside the normal and Chern-insulator phase 
depend slightly on the direction, they become universal 
at (A/t2)cr- Again, this is a signature of the metalhc 
character. 

It is interesting that the power of the power-law decay 
at the phase boundary seems to be exactly an integer and 
that it differs by 1 for different . This behavior can be 
understood in the following way: ^c«/3(R') of Eq. (23) is es- 
sentially the Fourier transform of the kernel and 
it is well known that discontinuities in the kernel deter- 
mine the decay behavior of the resulting quantity. In one 
dimension the discontinuities are related to the decay like 
^-C+i)^ where / is the number of continuous derivatives 
of the kernel. ^^'^^ Unfortunately, in two dimensions the 
situation is more complex and the resulting BZ integrals 
cannot easily be solved analytically. Nevertheless, we 
give heuristic arguments that a similar expression holds 
for higher dimensions. 

To this end, we solve for analytic expressions of C^^ by 
diagonalizing the Hamiltonian iJ (k) in Eq. (4) . In turn, 
we find analytic expressions for the kernel C^^C^^. Next, 
we switch to polar coordinates k = {kx, ky) ^ k = {k,(j)), 
replace A/t2 by (A/f2)cr = Sy/Ssimp, and expand the 
kernel around the Dirac point in orders of k: 

k,^k _ 1 3t2sin^ 

L/ii Oil — ■—■ ; K 

5t2 cos 3(j) sin ip 2 
~t~ ^ k 

+ 0{k^) (24) 
^k*^k _ (1 - zV3)e-'^ 

*-^ll<-^12 — ^ 

(3 -I- i\/3)e-''^ sin3^!> 
+ 48 " 

+ 0(fc2) (25) 

From the above expansions it is apparent that Cj'i'C}^ 
has its first discontinuity in first order in k. Hence, there 

are I = 1 continuous derivatives. On the other hand, due 
to the e^*"^ term, C\;iC^2 has already a discontinuity in 
zeroth order in fc, i.e. / = 0. This is consistent with 
the numerical results for ^a/3(R) in Fig. 9 if we assume 
that the decay in two dimensions is according to 



Equations (25) and (24) are thus consistent with a decay 
of R^^ and R^^ , respectively. 

In summary, the numerical and analytical arguments 
arc consistent in supporting the conclusion that the di- 
agonal and off-diagonal elements of ^q,^(R) decay as R~^ 
and R^^ respectively. An arbitrary pair of coordinates r 
and r' in Eq. (22) will involve a linear combination of con- 
tributions coming from diagonal and off-diagonal terms, 
so the final conclusion is that the decay of the density 
matrix will be as R~^ exactly on the NI/CI boundary, 
and exponential for any point lying within the normal- 
insulator or Chern-insulator phase. 

Above, we have evaluated the density matrix p(r, r') 
for periodic samples. For finite samples, we expect a 
parallel behavior to hold for points r and r' deep inside 
the bulk. However, if both points are chosen to be near 
the surface of a Chern-insulator sample, one may expect 
that the presence of metallic chiral edge states will induce 
a power-law decay with the distance between r and r' as 
measured along the perimeter. Preliminary calculations 
on finite samples appear consistent with this picture. 



VII. CONCLUSIONS 

Wc have performed numerical and analytical calcula- 
tions to study the behavior of several properties of the 
Haldane model as the system undergoes a transition from 
the normal-insulator phase to the Chern-insulator phase. 
We first showed how the usual methods of constructing 
Wannier functions break down for Chern insulators. We 
then investigated several quantities related to electron 
localization. We found that the total spread functional, 
which is finite in normal insulators, diverges in the case of 
a Chern insulator. However, when the spread functional 
is decomposed into its gauge-independent and gauge- 
dependent parts, the former is found to remain finite in 
a Chern insulator, while only the latter diverges. The lo- 
calization length of Resta and Sorella,* which is related 
to the gauge-independent part of the spread functional, 
thus remains finite. However, the localization length in- 
creases and diverges logarithmically as one approaches 
the NI/CI transition. Similarly, when inspecting the den- 
sity matrix, we find that it decays exponentially inside 
both the normal and Chern-insulator phases, but that 
the decay length increases as the phase boundary is ap- 
proached, and the behavior crosses over to a power-law 
decay exactly at the phase boundary. 

We thus find that a system that is sitting right on 
the NI/CI boundary has a kind of semimetallic charac- 
ter similar to that of graphene. in which the valence and 
conduction bands touch at one (for the Haldane model) 
or two (for graphene) Dirac points in the BZ. When the 
system is in the Chern-insulator phase, it still has rem- 
nants of metallic behavior in the presence of metallic edge 
states, the divergence of the total spread functional, and 
the difficulty of constructing Wannier functions. 

Our results were obtained here for a specific realiza- 
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tion of a Chern insulator, namely, the Haldane model. 
While it seems very likely that the localization properties 
found here will apply to other Chern-insulator systems, 
it remains to test this hypothesis by carrying out similar 
studies on other systems. It would also be of considerable 
interest to extend the current study to three-dimensional 
Chern-insulator crystals, and to continuum, as opposed 
to tight-binding, models. These could be fruitful avenues 
for future investigations. 
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